{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9c7526c4-3a64-4838-be44-86fb9c67dfb0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 3 folds for each of 324 candidates, totalling 972 fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Program Files\\python\\lib\\site-packages\\sklearn\\model_selection\\_validation.py:547: FitFailedWarning: \n",
      "324 fits failed out of a total of 972.\n",
      "The score on these train-test partitions for these parameters will be set to nan.\n",
      "If these failures are not expected, you can try to debug them by setting error_score='raise'.\n",
      "\n",
      "Below are more details about the failures:\n",
      "--------------------------------------------------------------------------------\n",
      "240 fits failed with the following error:\n",
      "Traceback (most recent call last):\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\model_selection\\_validation.py\", line 895, in _fit_and_score\n",
      "    estimator.fit(X_train, y_train, **fit_params)\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\base.py\", line 1467, in wrapper\n",
      "    estimator._validate_params()\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\base.py\", line 666, in _validate_params\n",
      "    validate_parameter_constraints(\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\utils\\_param_validation.py\", line 95, in validate_parameter_constraints\n",
      "    raise InvalidParameterError(\n",
      "sklearn.utils._param_validation.InvalidParameterError: The 'max_features' parameter of RandomForestRegressor must be an int in the range [1, inf), a float in the range (0.0, 1.0], a str among {'sqrt', 'log2'} or None. Got 'auto' instead.\n",
      "\n",
      "--------------------------------------------------------------------------------\n",
      "84 fits failed with the following error:\n",
      "Traceback (most recent call last):\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\model_selection\\_validation.py\", line 895, in _fit_and_score\n",
      "    estimator.fit(X_train, y_train, **fit_params)\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\base.py\", line 1467, in wrapper\n",
      "    estimator._validate_params()\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\base.py\", line 666, in _validate_params\n",
      "    validate_parameter_constraints(\n",
      "  File \"E:\\Program Files\\python\\lib\\site-packages\\sklearn\\utils\\_param_validation.py\", line 95, in validate_parameter_constraints\n",
      "    raise InvalidParameterError(\n",
      "sklearn.utils._param_validation.InvalidParameterError: The 'max_features' parameter of RandomForestRegressor must be an int in the range [1, inf), a float in the range (0.0, 1.0], a str among {'log2', 'sqrt'} or None. Got 'auto' instead.\n",
      "\n",
      "  warnings.warn(some_fits_failed_message, FitFailedWarning)\n",
      "E:\\Program Files\\python\\lib\\site-packages\\sklearn\\model_selection\\_search.py:1051: UserWarning: One or more of the test scores are non-finite: [       nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan 0.91119561 0.9096049  0.90846846\n",
      " 0.90594963 0.90171645 0.90215624 0.88889693 0.88542181 0.88542137\n",
      " 0.89986323 0.89870243 0.89980097 0.89968114 0.8970969  0.89765706\n",
      " 0.88450613 0.88259688 0.88330675 0.87689484 0.87773039 0.87753408\n",
      " 0.87689484 0.87773039 0.87753408 0.87603634 0.87242323 0.87333008\n",
      " 0.92279688 0.92041805 0.92082444 0.92375289 0.92081017 0.92133341\n",
      " 0.91801582 0.9148656  0.91567926 0.9255906  0.92198427 0.92207074\n",
      " 0.92139653 0.91996383 0.92084565 0.91286662 0.91234174 0.91382133\n",
      " 0.9069372  0.90585684 0.90633947 0.9069372  0.90585684 0.90633947\n",
      " 0.90332666 0.90338376 0.90468997        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      " 0.90816236 0.90800884 0.90700802 0.90578996 0.90163728 0.90209855\n",
      " 0.88889693 0.88542181 0.88542137 0.89987144 0.89870656 0.89980367\n",
      " 0.89968156 0.897097   0.89765708 0.88450613 0.88259688 0.88330675\n",
      " 0.87689484 0.87773039 0.87753408 0.87689484 0.87773039 0.87753408\n",
      " 0.87603634 0.87242323 0.87333008 0.92209813 0.91987614 0.92075716\n",
      " 0.92373472 0.92082844 0.92135048 0.91801582 0.9148656  0.91567926\n",
      " 0.9255906  0.92198427 0.92207074 0.92139653 0.91996383 0.92084565\n",
      " 0.91286662 0.91234174 0.91382133 0.9069372  0.90585684 0.90633947\n",
      " 0.9069372  0.90585684 0.90633947 0.90332666 0.90338376 0.90468997\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan 0.90816236 0.90800884 0.90700802\n",
      " 0.90578996 0.90163728 0.90209855 0.88889693 0.88542181 0.88542137\n",
      " 0.89987144 0.89870656 0.89980367 0.89968156 0.897097   0.89765708\n",
      " 0.88450613 0.88259688 0.88330675 0.87689484 0.87773039 0.87753408\n",
      " 0.87689484 0.87773039 0.87753408 0.87603634 0.87242323 0.87333008\n",
      " 0.92209813 0.91987614 0.92075716 0.92373472 0.92082844 0.92135048\n",
      " 0.91801582 0.9148656  0.91567926 0.9255906  0.92198427 0.92207074\n",
      " 0.92139653 0.91996383 0.92084565 0.91286662 0.91234174 0.91382133\n",
      " 0.9069372  0.90585684 0.90633947 0.9069372  0.90585684 0.90633947\n",
      " 0.90332666 0.90338376 0.90468997        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      "        nan        nan        nan        nan        nan        nan\n",
      " 0.90816236 0.90800884 0.90700802 0.90578996 0.90163728 0.90209855\n",
      " 0.88889693 0.88542181 0.88542137 0.89987144 0.89870656 0.89980367\n",
      " 0.89968156 0.897097   0.89765708 0.88450613 0.88259688 0.88330675\n",
      " 0.87689484 0.87773039 0.87753408 0.87689484 0.87773039 0.87753408\n",
      " 0.87603634 0.87242323 0.87333008 0.92209813 0.91987614 0.92075716\n",
      " 0.92373472 0.92082844 0.92135048 0.91801582 0.9148656  0.91567926\n",
      " 0.9255906  0.92198427 0.92207074 0.92139653 0.91996383 0.92084565\n",
      " 0.91286662 0.91234174 0.91382133 0.9069372  0.90585684 0.90633947\n",
      " 0.9069372  0.90585684 0.90633947 0.90332666 0.90338376 0.90468997]\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "🎯 最佳参数组合: {'max_depth': 10, 'max_features': 'log2', 'min_samples_leaf': 2, 'min_samples_split': 2, 'n_estimators': 100}\n",
      "⭐️ 最佳得分 R²: 0.9256\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import rasterio\n",
    "from rasterio.enums import Resampling\n",
    "from rasterio.warp import reproject\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import Normalize\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from skimage.measure import block_reduce\n",
    "from io import BytesIO\n",
    "import base64\n",
    "import os\n",
    "\n",
    "# 创建保存输出目录\n",
    "os.makedirs(\"output\", exist_ok=True)\n",
    "\n",
    "# 读取Landsat影像\n",
    "def read_raster(file_path):\n",
    "    with rasterio.open(file_path) as src:\n",
    "        data = src.read(1)\n",
    "        transform = src.transform\n",
    "        crs = src.crs\n",
    "        nodata = src.nodata\n",
    "    return data, transform, crs, nodata\n",
    "\n",
    "# 重采样\n",
    "def resample_raster(data, transform, target_shape):\n",
    "    new_transform = transform * transform.scale(\n",
    "        (data.shape[1] / target_shape[1]), (data.shape[0] / target_shape[0])\n",
    "    )\n",
    "    resampled_data = np.zeros(target_shape, dtype=np.float32)\n",
    "    reproject(\n",
    "        data, resampled_data,\n",
    "        src_transform=transform, src_crs='EPSG:4326',\n",
    "        dst_transform=new_transform, dst_crs='EPSG:4326',\n",
    "        resampling=Resampling.nearest\n",
    "    )\n",
    "    return resampled_data, new_transform\n",
    "\n",
    "# 降采样\n",
    "def aggregate_to_lower_resolution(data, factor=5, method='mean'):\n",
    "    if method == 'mean':\n",
    "        return block_reduce(data, block_size=(factor, factor), func=np.mean)\n",
    "    elif method == 'max':\n",
    "        return block_reduce(data, block_size=(factor, factor), func=np.max)\n",
    "    elif method == 'min':\n",
    "        return block_reduce(data, block_size=(factor, factor), func=np.min)\n",
    "    else:\n",
    "        raise ValueError(\"method must be 'mean', 'max', or 'min'\")\n",
    "\n",
    "# NDVI\n",
    "def calculate_ndvi(nir_band, red_band):\n",
    "    ndvi = (nir_band - red_band) / (nir_band + red_band + 1e-5)\n",
    "    return ndvi\n",
    "\n",
    "# 调参\n",
    "def train_model_with_tuning(features, target):\n",
    "    param_grid = {\n",
    "        'n_estimators': [100, 200, 300],\n",
    "        'max_depth': [10, 20, 30, None],\n",
    "        'min_samples_split': [2, 5, 10],\n",
    "        'min_samples_leaf': [1, 2, 4],\n",
    "        'max_features': ['auto', 'sqrt', 'log2']\n",
    "    }\n",
    "    model = RandomForestRegressor(random_state=42)\n",
    "    grid_search = GridSearchCV(model, param_grid, cv=3, scoring='r2', n_jobs=-1, verbose=2)\n",
    "    grid_search.fit(features, target)\n",
    "    print(\"🎯 最佳参数组合:\", grid_search.best_params_)\n",
    "    print(f\"⭐️ 最佳得分 R²: {grid_search.best_score_:.4f}\")\n",
    "    return grid_search.best_estimator_\n",
    "\n",
    "# 模型评估\n",
    "def evaluate_model(model, features, target):\n",
    "    predictions = model.predict(features)\n",
    "    mse = mean_squared_error(target, predictions)\n",
    "    r2 = r2_score(target, predictions)\n",
    "    return predictions, mse, r2\n",
    "\n",
    "# 保存图为 base64 字符串\n",
    "def save_plot_to_base64():\n",
    "    buf = BytesIO()\n",
    "    plt.savefig(buf, format=\"png\", bbox_inches='tight')\n",
    "    buf.seek(0)\n",
    "    base64_str = base64.b64encode(buf.read()).decode(\"utf-8\")\n",
    "    buf.close()\n",
    "    plt.close()\n",
    "    return base64_str\n",
    "\n",
    "# NDVI 图\n",
    "def generate_ndvi_plots(ndvi_100m, ndvi_500m):\n",
    "    img_list = []\n",
    "\n",
    "    plt.imshow(ndvi_100m, cmap='RdYlGn')\n",
    "    plt.title(\"NDVI (100m)\")\n",
    "    plt.colorbar()\n",
    "    img_list.append(save_plot_to_base64())\n",
    "\n",
    "    plt.imshow(ndvi_500m, cmap='RdYlGn')\n",
    "    plt.title(\"NDVI Aggregated to 500m\")\n",
    "    plt.colorbar()\n",
    "    img_list.append(save_plot_to_base64())\n",
    "\n",
    "    html = \"<html><body><h1>NDVI 可视化</h1>\"\n",
    "    for img in img_list:\n",
    "        html += f'<img src=\"data:image/png;base64,{img}\" style=\"width:90%;\"><br>'\n",
    "    html += \"</body></html>\"\n",
    "\n",
    "    with open(\"output/ndvi_plots.html\", \"w\", encoding='utf-8') as f:\n",
    "        f.write(html)\n",
    "\n",
    "# 各特征-预测值关系图\n",
    "def generate_feature_plots(features, predictions, r2, mse):\n",
    "    html = \"<html><body><h1>特征与预测产量关系图</h1>\"\n",
    "    for feature_name in features.columns:\n",
    "        # 选择每个特征的列作为 x\n",
    "        x = features[feature_name].values\n",
    "        # 确保 y 是预测结果，并且对应的长度匹配\n",
    "        y = predictions\n",
    "\n",
    "        # 计算回归线的斜率和截距\n",
    "        m, b = np.polyfit(x, y, 1)\n",
    "\n",
    "        plt.figure(figsize=(6, 4))\n",
    "        plt.scatter(x, y, alpha=0.5)\n",
    "        # 绘制回归线\n",
    "        plt.plot(x, m * x + b, color='red', linestyle='--', label=f\"y = {m:.2f}x + {b:.2f}\")\n",
    "\n",
    "        # 添加标题和坐标轴标签\n",
    "        plt.title(f\"{feature_name} vs Predicted Yield\")\n",
    "        plt.xlabel(feature_name)\n",
    "        plt.ylabel(\"Predicted Yield\")\n",
    "\n",
    "        # 在图表上添加回归线方程\n",
    "        plt.text(0.95, 0.95, f'y = {m:.2f}x + {b:.2f}', \n",
    "                 transform=plt.gca().transAxes,\n",
    "                 verticalalignment='top', horizontalalignment='right', fontsize=10, color='blue')\n",
    "\n",
    "        img = save_plot_to_base64()\n",
    "        html += f'<img src=\"data:image/png;base64,{img}\" style=\"width:85%;\"><br>'\n",
    "    html += \"</body></html>\"\n",
    "    \n",
    "    with open(\"output/feature_plots.html\", \"w\", encoding='utf-8') as f:\n",
    "        f.write(html)\n",
    "\n",
    "\n",
    "# 产量随时间变化图\n",
    "def generate_yield_time_plot(yield_data):\n",
    "    days = np.arange(1, len(yield_data) + 1)\n",
    "    plt.figure(figsize=(8, 5))\n",
    "    plt.plot(days, yield_data, marker='o', color='green')\n",
    "    plt.xlabel(\"Days\")\n",
    "    plt.ylabel(\"Predicted Yield (Kg/Ha)\")\n",
    "    plt.title(\"Predicted Corn Yield Over Time\")\n",
    "    plt.grid(True)\n",
    "\n",
    "    img = save_plot_to_base64()\n",
    "    html = f\"\"\"<html><body><h1>产量时间预测图</h1>\n",
    "               <img src=\"data:image/png;base64,{img}\" style=\"width:90%;\">\n",
    "               </body></html>\"\"\"\n",
    "    with open(\"output/yield_time_plot.html\", \"w\", encoding='utf-8') as f:\n",
    "        f.write(html)\n",
    "\n",
    "# 主流程\n",
    "def main():\n",
    "    # 读取影像\n",
    "    nir_data_100m, transform_100m, _, _ = read_raster(\"Landsat_guan_100m_2023.tif\")\n",
    "    red_data_500m, transform_500m, _, _ = read_raster(\"Landsat_guan_500m_2023.tif\")\n",
    "\n",
    "    # 重采样红波段\n",
    "    red_resampled, _ = resample_raster(red_data_500m, transform_500m, nir_data_100m.shape)\n",
    "\n",
    "    # NDVI计算与降采样\n",
    "    ndvi_100m = calculate_ndvi(nir_data_100m, red_resampled)\n",
    "    ndvi_500m = aggregate_to_lower_resolution(ndvi_100m, factor=5, method='mean')\n",
    "\n",
    "    # NDVI 图合成网页\n",
    "    generate_ndvi_plots(ndvi_100m, ndvi_500m)\n",
    "\n",
    "    # 加载数据\n",
    "    df = pd.read_csv(\"simulated_corn_yield.csv\")\n",
    "    features = df[['Temperature', 'Precipitation', 'Humidity', 'WindSpeed',\n",
    "                   'SoilMoisture', 'SoilpH', 'FlowRate', 'WaterLevel']]\n",
    "    target = df['Yield']\n",
    "\n",
    "    # 模型训练与预测\n",
    "    model = train_model_with_tuning(features, target)\n",
    "    predictions, mse, r2 = evaluate_model(model, features, target)\n",
    "\n",
    "    # 保存各特征图网页\n",
    "    generate_feature_plots(features, predictions, r2, mse)\n",
    "\n",
    "    # 预测产量时间图网页\n",
    "    generate_yield_time_plot(predictions)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b8d094a-8091-4b6f-830b-21f7f4cfa58c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "914b3b71-a52b-453a-8b49-3263e3cb5d5f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
